function res = CampbellSchillerRegression(Data,m)

% The CS regression reads:
% y_t+1_(k-m) - y_t_k = alpha + betta*m/(k-m)*(y_t_k-y_t_m) + eps_t+1
%
% INPUT: Must be a data matrix of dim: T * numYields, 
% where yields have maturities of 1 quarter, 2 quarters,3 quarters, ...
[T,n]          = size(Data);
matSelect      = m:m:n;
shorRate       = Data(1:T-m,m);
yieldsUsed     = Data(:,matSelect);
numYields      = size(yieldsUsed,2);
CSBetta        = NaN(2,numYields);
CSBetta_tstat  = NaN(2,numYields);
CSBetta_se     = NaN(2,numYields);
CSBetta_seBoot = NaN(2,numYields);
Ydata          = NaN(T-m,numYields);
Xdata          = NaN(T-m,numYields);
R2             = NaN(1,numYields);
for i=1:numYields-1
   Y  = yieldsUsed(1+m:T,i)-yieldsUsed(1:T-m,i+1);
   %x  = (yieldsUsed(1:T-m,i+1)-shorRate)*m/(matSelect(1,i+1)-m);
   x  = (yieldsUsed(1:T-m,i+1)-shorRate);
   X  = [ones(T-m,1) x*m/(matSelect(1,i+1)-m)];
   resOLS            = nwest(Y,X,m+1);
   CSBetta(:,i+1)    = resOLS.beta;
   CSBetta_se(:,i+1) = resOLS.beta./resOLS.tstat;
   CSBetta_tstat(:,i+1) = resOLS.tstat;
   Ydata(:,i+1)      = Y;
   Xdata(:,i+1)      = x*m/(matSelect(1,i+1)-m);
   R2(:,i+1)         = resOLS.rsqr;
   % Bootstrapping the SE
   if T < 1000
       bootWidth = 10;
       numBoots  = 10000;
       %numBoots  = 1D5;
       rng(1,'twister')
       bsdata=block_bootstrap([Y,X],bootWidth,numBoots);
       betaBoot  = zeros(2,numBoots);
       for s=1:numBoots
           dataBoot = bsdata(:,:,s);
           Yboot    = dataBoot(:,1);
           Xboot    = dataBoot(:,2:3);
           betaBoot(:,s) = (Xboot'*Xboot)\(Xboot'*Yboot);
       end
       CSBetta_seBoot(:,i+1) = std(betaBoot,0,2);
   end
end

% Output
res.CSbetta        = CSBetta;
res.maturities     = matSelect;
res.CSbetta_se     = CSBetta_se;
res.CSBetta_tstat  = CSBetta_tstat;
res.CSBetta_CI95   = [CSBetta(2,:)-1.96*CSBetta_se(2,:);CSBetta(2,:)+1.96*CSBetta_se(2,:)];
res.CSBetta_CI99   = [CSBetta(2,:)-2.575*CSBetta_se(2,:);CSBetta(2,:)+2.575*CSBetta_se(2,:)];
if T < 1000
    res.CSBetta_seBoot = CSBetta_seBoot;
    res.CSBetta_CI95Boot =[CSBetta(2,:)-1.96*CSBetta_seBoot(2,:);CSBetta(2,:)+1.96*CSBetta_seBoot(2,:)];
end
res.Ydata          = Ydata;
res.Xdata          = Xdata;
res.R2             = R2;


